function p = gaussian_pdf( sigma, Nsyn, neurons_in_reach, max_r  )

% determine the probability of forming a connection.

total = 0;
pos = [];
gsum = 0;

scale = Nsyn*(2*max_r)/myerf(sigma,max_r)/neurons_in_reach
uniform = false;
if(scale > 1/gaussian_1d(sigma,0))
    disp('warning: not enough cells in reach to get gaussian pdf; switching to uniform distribution');
    uniform = true;
end
for i = 1:neurons_in_reach

    x=2*max_r*rand-max_r;

    if uniform
        g= Nsyn/neurons_in_reach;
    else
        g= scale * gaussian_1d(sigma, x);
    end
    
    
    gsum = gsum + g;
    
    s = rand < g;
   
    
    total = total + s;

    if s
        pos = [pos x];
    end
end

disp( ['Wanted ' num2str(Nsyn) '; got ' num2str(total)]);
disp( [ 'gsum was ' num2str(gsum) ]);

hist(pos);




function xa = gaussian_1d(sigma_p,x);

x0 = 0;


xa = 1/(sigma_p * sqrt(2*pi)) * exp ((-(-x-x0)^2) / (2*sigma_p^2));
